home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / DJINC106.ARJ / MATH-688.H < prev    next >
C/C++ Source or Header  |  1992-03-29  |  8KB  |  465 lines
  1. /******************************************************************\
  2. *                                   *
  3. *  <math-68881.h>        last modified: 18 May 1989.       *
  4. *                                   *
  5. *  Copyright (C) 1989 by Matthew Self.                   *
  6. *  You may freely distribute verbatim copies of this software       *
  7. *  provided that this copyright notice is retained in all copies.  *
  8. *  You may distribute modifications to this software under the     *
  9. *  conditions above if you also clearly note such modifications    *
  10. *  with their author and date.                               *
  11. *                                   *
  12. *  Note:  errno is not set to EDOM when domain errors occur for    *
  13. *  most of these functions.  Rather, it is assumed that the       *
  14. *  68881's OPERR exception will be enabled and handled           *
  15. *  appropriately by the    operating system.  Similarly, overflow       *
  16. *  and underflow do not set errno to ERANGE.               *
  17. *                                   *
  18. *  Send bugs to Matthew Self (self@bayes.arc.nasa.gov).           *
  19. *                                   *
  20. \******************************************************************/
  21.  
  22. #include <errno.h>
  23.  
  24. #ifndef HUGE_VAL
  25. #define HUGE_VAL                            \
  26. ({                                    \
  27.   double huge_val;                            \
  28.                                     \
  29.   __asm ("fmove%.d #0x7ff0000000000000,%0"    /* Infinity */        \
  30.      : "=f" (huge_val)                        \
  31.      : /* no inputs */);                        \
  32.   huge_val;                                \
  33. })
  34. #endif
  35.  
  36. __inline static const double sin (double x)
  37. {
  38.   double value;
  39.  
  40.   __asm ("fsin%.x %1,%0"
  41.      : "=f" (value)
  42.      : "f" (x));
  43.   return value;
  44. }
  45.  
  46. __inline static const double cos (double x)
  47. {
  48.   double value;
  49.  
  50.   __asm ("fcos%.x %1,%0"
  51.      : "=f" (value)
  52.      : "f" (x));
  53.   return value;
  54. }
  55.  
  56. __inline static const double tan (double x)
  57. {
  58.   double value;
  59.  
  60.   __asm ("ftan%.x %1,%0"
  61.      : "=f" (value)
  62.      : "f" (x));
  63.   return value;
  64. }
  65.  
  66. __inline static const double asin (double x)
  67. {
  68.   double value;
  69.  
  70.   __asm ("fasin%.x %1,%0"
  71.      : "=f" (value)
  72.      : "f" (x));
  73.   return value;
  74. }
  75.  
  76. __inline static const double acos (double x)
  77. {
  78.   double value;
  79.  
  80.   __asm ("facos%.x %1,%0"
  81.      : "=f" (value)
  82.      : "f" (x));
  83.   return value;
  84. }
  85.  
  86. __inline static const double atan (double x)
  87. {
  88.   double value;
  89.  
  90.   __asm ("fatan%.x %1,%0"
  91.      : "=f" (value)
  92.      : "f" (x));
  93.   return value;
  94. }
  95.  
  96. __inline static const double atan2 (double y, double x)
  97. {
  98.   double pi, pi_over_2;
  99.  
  100.   __asm ("fmovecr%.x %#0,%0"        /* extended precision pi */
  101.      : "=f" (pi)
  102.      : /* no inputs */ );
  103.   __asm ("fscale%.b %#-1,%0"        /* no loss of accuracy */
  104.      : "=f" (pi_over_2)
  105.      : "0" (pi));
  106.   if (x > 0)
  107.     {
  108.       if (y > 0)
  109.     {
  110.       if (x > y)
  111.         return atan (y / x);
  112.       else
  113.         return pi_over_2 - atan (x / y);
  114.     }
  115.       else
  116.     {
  117.       if (x > -y)
  118.         return atan (y / x);
  119.       else
  120.         return - pi_over_2 - atan (x / y);
  121.     }
  122.     }
  123.   else
  124.     {
  125.       if (y > 0)
  126.     {
  127.       if (-x > y)
  128.         return pi + atan (y / x);
  129.       else
  130.         return pi_over_2 - atan (x / y);
  131.     }
  132.       else
  133.     {
  134.       if (-x > -y)
  135.         return - pi + atan (y / x);
  136.       else if (y < 0)
  137.         return - pi_over_2 - atan (x / y);
  138.       else
  139.         {
  140.           double value;
  141.  
  142.           errno = EDOM;
  143.           __asm ("fmove%.d %#0rnan,%0"     /* quiet NaN */
  144.              : "=f" (value)
  145.              : /* no inputs */);
  146.           return value;
  147.         }
  148.     }
  149.     }
  150. }
  151.  
  152. __inline static const double sinh (double x)
  153. {
  154.   double value;
  155.  
  156.   __asm ("fsinh%.x %1,%0"
  157.      : "=f" (value)
  158.      : "f" (x));
  159.   return value;
  160. }
  161.  
  162. __inline static const double cosh (double x)
  163. {
  164.   double value;
  165.  
  166.   __asm ("fcosh%.x %1,%0"
  167.      : "=f" (value)
  168.      : "f" (x));
  169.   return value;
  170. }
  171.  
  172. __inline static const double tanh (double x)
  173. {
  174.   double value;
  175.  
  176.   __asm ("ftanh%.x %1,%0"
  177.      : "=f" (value)
  178.      : "f" (x));
  179.   return value;
  180. }
  181.  
  182. __inline static const double atanh (double x)
  183. {
  184.   double value;
  185.  
  186.   __asm ("fatanh%.x %1,%0"
  187.      : "=f" (value)
  188.      : "f" (x));
  189.   return value;
  190. }
  191.  
  192. __inline static const double exp (double x)
  193. {
  194.   double value;
  195.  
  196.   __asm ("fetox%.x %1,%0"
  197.      : "=f" (value)
  198.      : "f" (x));
  199.   return value;
  200. }
  201.  
  202. __inline static const double expm1 (double x)
  203. {
  204.   double value;
  205.  
  206.   __asm ("fetoxm1%.x %1,%0"
  207.      : "=f" (value)
  208.      : "f" (x));
  209.   return value;
  210. }
  211.  
  212. __inline static const double log (double x)
  213. {
  214.   double value;
  215.  
  216.   __asm ("flogn%.x %1,%0"
  217.      : "=f" (value)
  218.      : "f" (x));
  219.   return value;
  220. }
  221.  
  222. __inline static const double log1p (double x)
  223. {
  224.   double value;
  225.  
  226.   __asm ("flognp1%.x %1,%0"
  227.      : "=f" (value)
  228.      : "f" (x));
  229.   return value;
  230. }
  231.  
  232. __inline static const double log10 (double x)
  233. {
  234.   double value;
  235.  
  236.   __asm ("flog10%.x %1,%0"
  237.      : "=f" (value)
  238.      : "f" (x));
  239.   return value;
  240. }
  241.  
  242. __inline static const double sqrt (double x)
  243. {
  244.   double value;
  245.  
  246.   __asm ("fsqrt%.x %1,%0"
  247.      : "=f" (value)
  248.      : "f" (x));
  249.   return value;
  250. }
  251.  
  252. __inline static const double pow (const double x, const double y)
  253. {
  254.   if (x > 0)
  255.     return exp (y * log (x));
  256.   else if (x == 0)
  257.     {
  258.       if (y > 0)
  259.     return 0.0;
  260.       else
  261.     {
  262.       double value;
  263.  
  264.       errno = EDOM;
  265.       __asm ("fmove%.d %#0rnan,%0"        /* quiet NaN */
  266.          : "=f" (value)
  267.          : /* no inputs */);
  268.       return value;
  269.     }
  270.     }
  271.   else
  272.     {
  273.       double temp;
  274.  
  275.       __asm ("fintrz%.x %1,%0"
  276.          : "=f" (temp)            /* integer-valued float */
  277.          : "f" (y));
  278.       if (y == temp)
  279.         {
  280.       int i = (int) y;
  281.       
  282.       if ((i & 1) == 0)            /* even */
  283.         return exp (y * log (x));
  284.       else
  285.         return - exp (y * log (x));
  286.         }
  287.       else
  288.         {
  289.       double value;
  290.  
  291.       errno = EDOM;
  292.       __asm ("fmove%.d %#0rnan,%0"        /* quiet NaN */
  293.          : "=f" (value)
  294.          : /* no inputs */);
  295.       return value;
  296.         }
  297.     }
  298. }
  299.  
  300. __inline static const double fabs (double x)
  301. {
  302.   double value;
  303.  
  304.   __asm ("fabs%.x %1,%0"
  305.      : "=f" (value)
  306.      : "f" (x));
  307.   return value;
  308. }
  309.  
  310. __inline static const double ceil (double x)
  311. {
  312.   int rounding_mode, round_up;
  313.   double value;
  314.  
  315.   __asm volatile ("fmove%.l fpcr,%0"
  316.           : "=dm" (rounding_mode)
  317.           : /* no inputs */ );
  318.   round_up = rounding_mode | 0x30;
  319.   __asm volatile ("fmove%.l %0,fpcr"
  320.           : /* no outputs */
  321.           : "dmi" (round_up));
  322.   __asm volatile ("fint%.x %1,%0"
  323.           : "=f" (value)
  324.           : "f" (x));
  325.   __asm volatile ("fmove%.l %0,fpcr"
  326.           : /* no outputs */
  327.           : "dmi" (rounding_mode));
  328.   return value;
  329. }
  330.  
  331. __inline static const double floor (double x)
  332. {
  333.   int rounding_mode, round_down;
  334.   double value;
  335.  
  336.   __asm volatile ("fmove%.l fpcr,%0"
  337.           : "=dm" (rounding_mode)
  338.           : /* no inputs */ );
  339.   round_down = (rounding_mode & ~0x10)
  340.         | 0x20;
  341.   __asm volatile ("fmove%.l %0,fpcr"
  342.           : /* no outputs */
  343.           : "dmi" (round_down));
  344.   __asm volatile ("fint%.x %1,%0"
  345.           : "=f" (value)
  346.           : "f" (x));
  347.   __asm volatile ("fmove%.l %0,fpcr"
  348.           : /* no outputs */
  349.           : "dmi" (rounding_mode));
  350.   return value;
  351. }
  352.  
  353. __inline static const double rint (double x)
  354. {
  355.   int rounding_mode, round_nearest;
  356.   double value;
  357.  
  358.   __asm volatile ("fmove%.l fpcr,%0"
  359.           : "=dm" (rounding_mode)
  360.           : /* no inputs */ );
  361.   round_nearest = rounding_mode & ~0x30;
  362.   __asm volatile ("fmove%.l %0,fpcr"
  363.           : /* no outputs */
  364.           : "dmi" (round_nearest));
  365.   __asm volatile ("fint%.x %1,%0"
  366.           : "=f" (value)
  367.           : "f" (x));
  368.   __asm volatile ("fmove%.l %0,fpcr"
  369.           : /* no outputs */
  370.           : "dmi" (rounding_mode));
  371.   return value;
  372. }
  373.  
  374. __inline static const double fmod (double x, double y)
  375. {
  376.   double value;
  377.  
  378.   __asm ("fmod%.x %2,%0"
  379.      : "=f" (value)
  380.      : "0" (x),
  381.        "f" (y));
  382.   return value;
  383. }
  384.  
  385. __inline static const double drem (double x, double y)
  386. {
  387.   double value;
  388.  
  389.   __asm ("frem%.x %2,%0"
  390.      : "=f" (value)
  391.      : "0" (x),
  392.        "f" (y));
  393.   return value;
  394. }
  395.  
  396. __inline static const double scalb (double x, int n)
  397. {
  398.   double value;
  399.  
  400.   __asm ("fscale%.l %2,%0"
  401.      : "=f" (value)
  402.      : "0" (x),
  403.        "dmi" (n));
  404.   return value;
  405. }
  406.  
  407. __inline static double logb (double x)
  408. {
  409.   double exponent;
  410.  
  411.   __asm ("fgetexp%.x %1,%0"
  412.      : "=f" (exponent)
  413.      : "f" (x));
  414.   return exponent;
  415. }
  416.  
  417. __inline static const double ldexp (double x, int n)
  418. {
  419.   double value;
  420.  
  421.   __asm ("fscale%.l %2,%0"
  422.      : "=f" (value)
  423.      : "0" (x),
  424.        "dmi" (n));
  425.   return value;
  426. }
  427.  
  428. __inline static double frexp (double x, int *exp)
  429. {
  430.   double float_exponent;
  431.   int int_exponent;
  432.   double mantissa;
  433.  
  434.   __asm ("fgetexp%.x %1,%0"
  435.      : "=f" (float_exponent)     /* integer-valued float */
  436.      : "f" (x));
  437.   int_exponent = (int) float_exponent;
  438.   __asm ("fgetman%.x %1,%0"
  439.      : "=f" (mantissa)        /* 1.0 <= mantissa < 2.0 */
  440.      : "f" (x));
  441.   if (mantissa != 0)
  442.     {
  443.       __asm ("fscale%.b %#-1,%0"
  444.          : "=f" (mantissa)        /* mantissa /= 2.0 */
  445.          : "0" (mantissa));
  446.       int_exponent += 1;
  447.     }
  448.   *exp = int_exponent;
  449.   return mantissa;
  450. }
  451.  
  452. __inline static double modf (double x, double *ip)
  453. {
  454.   double temp;
  455.  
  456.   __asm ("fintrz%.x %1,%0"
  457.      : "=f" (temp)            /* integer-valued float */
  458.      : "f" (x));
  459.   *ip = temp;
  460.   return x - temp;
  461. }
  462.  
  463.  
  464.  
  465.